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Transport phenomena are studied for a binary (AB) alloy on a rigid square lattice with 
nearest-neighbor attraction between unlike particles, assuming a small concentration c v of 
vacancies V being present, to which A(B) particles can jump with rates Ta(^b) in the case 
where the nearest neighbor attractive energy €ab is negligible in comparison with the thermal 
energy kgT in the system. This model exhibits a continuous order-disorder transition for 
concentrations ca,cb = 1 — ca — cy in the range c^' < ca < c C A2> w ^h c Xi* = 0- ~ m * ~ 
cy)/2, — (1 + m * ~ c v)/2, m* sa 0.25, the maximum critical temperature occurring 
for c* = ca = cb = (1 — cy)/2, i.e. m* = 0. This phase transition belongs to the d = 2 
Ising universality class, demonstrated by a finite size scaling analysis. From a study of 
mean-square displacements of tagged particles, self-diffusion coefficients are deduced, while 
applying chemical potential gradients allow the estimation of Onsager coefficients. Analyzing 
finally the decay with time of sinusoidal concentration variations that were prepared as 
initial condition, also the interdiffusion coefficient is obtained as function of concentration 
and temperature. As in the random alloy case (i.e., a noninteracting ABV- model) no simple 
relation between self-diffusion and interdiffusion is found. Unlike this model mean field 
theory cannot describe interdiffusion, however, even if the necessary Onsager coefficients are 
estimated via simulation. 



I. INTRODUCTION 



Understanding of atomic transport in multicomponent solids has been a longstanding challenge 
. In particular, the problem of interdiffusion in binary metallic alloys (as well as other types 
of mixed crystals) is very intricate: there is a delicate interplay between kinetic aspects that have 
a complicated energetics (such as jump rates of the various kinds of atoms to available vacant 
sites) and effects due to non-random arrangement of these atoms on the lattice sites (a problem 



which needs to be considered in the framework of statistical thermodynamics [14L llRl. Even 
the simplistic limiting case of perfectly random occupation of the sites of a rigid perfect lattice 
by two atomic species (A,B) and a small fraction of vacancies (V), where one assumes constant 
jump rates Ta,Tb of the two types of atoms to the vacant sites (i.e. jump rates that do not 



depend on the occupation of the sites surrounding the vacant sites), is highly nontrivial 
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finds that neither the self-diffusion coefficients Da,Db nor the interdiffusion coefficient Di n t can 
be analytically reliably predicted, given Ta,^b and the average concentration ca,cb; nor does a 



simple relation between Da,Db and Di n t exist 
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Recently, attention has been focused on this problem because of several fascinating develop- 
ments: (i) Progress with the electronic structure calculations of vacancy formation energies, jump 
rates etc. as well as better understanding of short range order parameters in alloys puts the 
"first-principles" calculation of interdiffusion and self-diffusion coefficients in ordered solid alloys 
such as Al(i_ x }Li x within reach [13] . (ii) Progress with the atom-tracking scanning tunnelling 
microscopy observation of atomic motions in two-dimensional surface alloys such as In atoms in 
Cu(OOl) surfaces 10] or Pd atoms in Cu(001) surfaces [ll| has provided compelling direct evidence 
for the operation of vacancy mediated surface diffusion. This is a nontrivial result, since competing 
mechanisms (surface atoms leave the topmost atomic layer to become adatoms on top of this layer 
[isl |. or direct exchange between neighboring surface atoms, "assisted" by the free space above 
the topmost monolayer of atoms at the crystal surface) can not be ruled out a priori. Of course, 
this finding enforces the hypothesis that vacancy mechanism dominates self- and interdiffusion 
processes in crystal lattices in the bulk 

In the present work we try to contribute to this problem, emphasizing the statistical mechanics 
approach by considering again a rigid lattice model but allowing for interactions causing a nontrivial 
long range order (or, at least, short range order) between the atoms in the system. We are not 
addressing a specific material, but rather try to elucidate the generic phenomena caused by the 
interplay of local correlations in the occupancy of lattice sites and the disparities in the jump rates 
r,4 and r# of the two species. Thus, our model is close in spirit to the work in Ref. [ljj and 
employs a related Monte Carlo simulation methodology Unlike ^?j], the present model does 
include a nearest-neighbor attraction between unlike neighbors, and thus nontrivial static order- 
disorder phenomena occur. As expected, we shall demonstrate that the resulting correlations 
in the occupancy of the lattice sites have a drastic effect on the transport phenomena, and hence 
cannot be neglected when one tries to interpret real data. We also emphasize that these correlation 
phenomena need a treatment beyond mean field level. We point out this fact, because sometimes a 
first principle electronic structure calculation is combined with statistical mechanics of mean field 
type or the cluster variation method 16], and such approximations then clearly spoil the desirable 
rigor. We also note that similar models as studied here have been frequently used to study domain 
growth in alloys that are quenched from the high temperature phase to a temperature below the 
order-disorder transition temperature j^ . 
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In Section II we describe our (two-dimensional) model. We restrict the present work to two- 
dimensional systems, since recently there has been great interest in two-dimensional alloys j^jj, and 
we hope that extensions of our modelling can make contact with corresponding experiments. In 
Section III we summarize our simulation methodology, while Section IV briefly reviews some per- 
tinent theoretical concepts and approximations we wish to test. Section V describes our numerical 
results, while Section VI summarizes our conclusions, and gives an outlook to future work. 



II. THE MODEL AND ITS STATIC PROPERTIES 



Having in mind the application of our work to two-dimensional surface alloys, we assume a 
perfect square lattice of adsorption sites (Fig. These adsorption sites can either be taken by 
an A-atom, a B-atom, or a vacancy. Therefore this model traditionally is also referred to as the 



ABV model 



ia 



221 ] . It can also be viewed as an extension of simple lattice gas models, where 



diffusion of a single species (A) occurs by hopping to vacant sites, to two components. Diffusion in 



lattice gases with a single species has been extensively studied 



p,yy, 
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but diffusion in a two-component lattice gas so far has been thoroughly examined only in the 
noninteracting case [l^]. Here we restrict attention to a model with strictly pairwise interactions 
between nearest neighbors only, which we denote as 6aa,£AB and ebb pairs. However, in general 
one can consider also energy parameters between pairs of lattice sites involving one (cav, cbv) or 
two (evv) vacancies, but here we do not consider the ABV model in full generality, but only the 
special case €av = £bv = Wv = 0, although from first principles electronic structure calculations 

n 

there is evidence that nonzero eav , £ bv > e V V ma y occur [22| . While all these parameters affect 
the diffusion behavior of the model, actually only a subset of them controls the static behavior. 
With respect to static properties of this model, the well-known transcription to the spin— 1 Blume- 
Emery-Griffiths model shows (see e.g. jsjj), that for constant concentrations only three interaction 
parameters would be needed. Note that although there are three concentration variables, ca,cb, 
and cy, due to the constraint ca + cb + cy = 1 only two of them are independent. Actually, the 
physically most interesting case is the limit cy — > 0, since in thermal equilibrium the concentration 
of vacancies is very small. In the noninteracting case |l7j |. it was found that many aspects of this 
limiting behavior cy — ► are already reproduced if the vacancy concentration is of the order of a 
few percent only, e.g. cy = 0.04, and in fact we shall adopt this choice in the case of the present 
simulations. Also for the interacting case the limit cy — > greatly simplifies matters, since then, 
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with respect to static properties, we have to consider only a single energy parameter e, defined by 



e = e A B ~ {e A A + ess)/2. (1) 

If e < 0, the model in thermal equilibrium will exhibit ordering, while for e > 0, phase separation 
occurs Q, 12 • I n the case °f a square lattice, the model in the limit cy — > is equivalent to 
the two-dimensional Ising model, for which some static properties of interest are exactly known 
In particular, for ca = cb = 1/2 the critical temperature T c is known exactly, namely 

k B T™ ax /\e\ = [ln(l + V2)}- 1 « 1.1345. (2) 

This is the maximum value of the critical temperature curve T c (ca) at which the order-disorder 
phase transition occurs. According to the well-known Bragg- Williams mean-field approximation, 
one would rather obtain ksT^ F /e = 2 than the result implied by Eq. (j2j), ksTc/e ~ 1.1345 
(37]]. Here and in the following, the maximum value T c (ca = 0.5) of the pure model without 
vacancies is simply denoted as T c . However, an even more important failure of the mean field 
theory is the prediction that an order-disorder transition from the disordered phase to a phase 
with long-range checkerboard order occurs over the entire concentration range, with T c (ca — > 0) — > 
0, T c (ca — ► 1) — ► 0, see jsTJ for a more detailed discussion of mean field theory. As a matter of 
'act, long range order is only possible for a much more restricted range of concentrations, namely 
3^ | 0.375 < ca < 0.625 (note that the pairwise character of the interactions implies a symmetry 



of the phase diagram around the line ca = 1/2, in the limit cy — > 14l. Ila. l33|). 

If we work with a small but nonzero concentration of vacancies cy, the maximum critical 
temperature no longer occurs at ca = cb = 1/2, but rather at c* = c A = eg = (1 — cy)/2, 
and the phase diagram is in this case symmetric around this concentration c*. Apart from this 
statement, there are no longer any exact results available, but it is fairly straightforward to obtain 
the phase diagram from standard Monte Carlo methods jly] with an accuracy that is sufficient for 
our purposes. Fig|5] shows our estimates of the phase boundary for cy = 0.04, in comparison with 



19, 



33 



3, 



such phase 



previous results for cy = 0. As has been well documented in the literature 
diagrams are conveniently mapped out by transforming the model to a magnetic Ising spin model 
(representing the cases that lattice site i is taken by an A— atom by spin up, B— atom by spin 
down, respectively) and considering the transition from the paramagnetic to the antiferromagnetic 
phase for various magnetic fields H (2 H = 11 a — fJ-B, if ^AA = £bb, and with \xa and \xb being 
the chemical potentials of A and B— particles, respectively). Estimating then the magnetization 
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m c (H) = m(H,T = T C (H)) at the phase boundary, one then obtains the corresponding critical 
concentrations from 

cT\T) = {l±m c {T)-c v )/2 (3) 

Fig. [21 shows that for cy = 0.04 the maximum critical temperature occurs for T c (cy = 0.04)/T c ~ 
0.905, and for T — > the phase boundary ends at the concentrations cj 4 * » 0.375, ~ 0.585. 
As it should be, the phase diagram is symmetric around c < ^ t ^ iax = (1 — cy)/2 = 0.48. The analysis 
indicates that the order-disorder transition stays second order throughout, also in the presence of 
this small vacancy concentration. Although it is clear that a vacancy concentration of cy = 0.04 
does have some clearly visible effects, in comparison to the model with cy — > 0, these changes do 
not affect the qualitative character of the phase behavior, but cause only minor modifications of 
quantitative details. For obtaining accurate results on the dynamic behavior of the model with 
a modest amount of computing time, working with sufficiently many vacancies on the lattice is 
mandatory. Note that for the diffusion studies we use a lattice of linear dimensions L = 1024, 
while the static phase diagram was extracted from a standard finite size scaling analysis |l9j], see 
Figs. 01 \M f° r an example, using sizes 24 < L < 192. Periodic boundary conditions are applied 
throughout. The static quantities that were analyzed in order to obtain the phase boundary are 
the antiferromagnetic order parameter ^ (we refer here to the transformation of the model to the 
Ising spin representation again) 

* = <M>, = L- 2 ^^(-l) fc+ ^ M , (4) 

k=l l=\ 

where k, I label the lattice sites in x— and y— direction, respectively. Similarly, the magnetization 
m is given by averaging all the spins without a phase factor 

L L 

m=(M), M = L- 2 Y J Y, S ^i ( 5 ) 

k=l l=\ 

and the susceptibility x an d staggered susceptibility x are obtained from the standard fluctuation 
relations 

X = L 2 ((M 2 ) - (M) 2 )/k B T , (6) 



X = L 2 {{<p 2 )-{\<p\) 2 )/k B T 



(7) 
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Note that in a finite system in the absence of symetry-breaking fields one needs to work with the 
average of the absolute value {\(ft\) rather than {(ft) in order to have a meaningful order parameter 



A further quantity useful for finding the location of the transition is the fourth order cumulant 
of the order parameter 3] 

U L = 1- {cft 4 )/[3{(ft 2 ) 2 } , (8) 

since the critical temperature can be found from the intersection of the cumulants plotted ver- 
sus temperature for different lattice sizes. For the two-dimensional Ising universality class, this 
intersection should occur for a value U* » 0.6107 |3j. 

Fig. |31 shows that this expectation is only rather roughly fulfilled. To some extent this may 
be attributed to statistical errors, but in addition probably for cy ^ there are somewhat larger 
corrections to finite size scaling than for the "pure" model (i.e., the model without vacancies). We 
have hence estimated T c {cy = 0.04) alternatively from a plot of the temperatures T C (L), where 
the maximum of x(T, L) for finite L occurs, versus the finite size scaling variable L~ l / U = L^ 1 
(remember u = 1 in the two-dimensional Ising model [si]]), see Fig. 0^. The quality of the finite 
size scaling "data collapse" of the order parameter (fig. IIJd) gives us confidence in the reliability of 
our procedures. 

We emphasize that the present paper concerns only the choice of the symmetric case, eaa = £bb- 
While any asymmetry between A and B, leading to €aa ^ £bb, has little effect on static properties 
for small ey, the distribution of the vacancies and their dynamics may get strongly affected by 



such an asymmetry 



20 



Finally, we mention a static quantity that plays a role in discussing the self-diffusion coefficient 



of particles in lattice gas models, the so-called "vacancy availability factor" 



flfl 



30] 



V = c v {l-a l ) (9) 



Here a\ is the standard Cowley- Warren short range order parameter ]14 , 15, 16, 331 for the nearest 
neighbor shell of a particle: a\ = if there is a random occupation of the lattice sites by any 
particles and vacancies (note that here we are not concerned with short range order describing the 
non-random occupation of A versus B particles on the lattice. Due to the symmetry €aa = £bb, 
there is also no need to consider separate vacancy availability factors for A and B particles). 
Actually, we expect that in the limit cy — > also ol\ — > 0, and then V = cy. Hence a calculation 
of q,\ can serve as a test whether the chosen vacancy concentration is small enough in order to 
reproduce the desired limit cy — > 0. 
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III. SIMULATION METHODOLOGY TO STUDY TRANSPORT PHENOMENA 

The Monte Carlo simulations consist of an initial part, necessary to equilibrate the system for 
the desired conditions, and a final part, where the transport coefficients of interest are "measured" 

n, 

in the simulation. While in the case of the completely random ABV alloy studied in [13] the 
generation of an initial configuration is straightforward, this is not so here, because depending on 
where the chosen state point (T, ca) is in the phase diagram, Fig. [2 we have long range order or not. 
If the system in equilibrium is in a state where long range order occurs, it is important to prepare 
the system in a monodomain sample: otherwise the presence of antiphase domain boundaries (3^] 
might spoil the results. In particular, at very low temperature interdiffusion could be strongly 
enhanced near such boundaries, in comparison with the bulk. Although such effects are interesting 
in their own right, they need separate study from bulk behavior, and are out of consideration here. 

Actually the best way to prepare the equilibrated initial configurations, in cases where long range 
order is present, is the use of the "magnetic" representation of the model as an Ising antiferromagnet 
in a field H (remember that H physically corresponds to the chemical potential difference between A 
and B particles Recording the magnetization m(T, H) as function of the field, one can choose 

the field such that states with the desired value of m and hence ca = (1 — cy — m)/2 result. The 
initial spin configuration is that of a perfect antiferromagnetic structure, from which a fraction cy 
of sites chosen at random is removed. The Monte Carlo algorithm that was used is the standard 
single spin flip Metropolis algorithm mixed with random exchanges of the vacancies with 
randomly chosen neighbors. Note that during this equilibration part the concentration ca is not 
strictly constant, but slightly fluctuating: this lack of conservation of ca is desirable, however, since 
"hydrodynamic slowing down" ^| of long wavelength concentration fluctuations would otherwise 
hamper the equilibration of concentration correlations at long distances. 

In the final stage of the Monte Carlo runs, of course, the spin flip Monte Carlo moves are shut 
off, since for the analysis of the diffusion constants the concentrations ca,cb = 1 — cy — ca need to 
be strictly conserved. Most straightforward is the estimation of the self-diffusion coefficients (also 
called '"tracer diffusion coefficients") Dt of tagged A and B particles, since there one simply can 
apply the Einstein relation 

(r 2 ) = 2dD t t, t^oo, r = fi(t)-fi{0), (10) 

d being the dimensionality of the lattice (d = 2 here), and fi{t) being the position of the i- 
th particle at time t. Fig. |5] illustrates the application of this method for a typical example, 
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in the case Ta/^b = 0.01, temperature T = 1.2 (in units of T c , Eq. |2J), and concentrations 
c A = 0.40, c B = 0.56, respectively. While the plot of (r 2 ) vs. t, for a total time of t = 10 4 MCS, 
look at first sight almost linear (Fig. left part), a closer look reveals a slight but systematic 
decrease of the slope of the (r 2 ) vs. t curve with increasing time. A similar observation was 
already reported by Kehr et al. 171 ]. who attributed this decrease of slope to the presence of a 
logarithmic correction. 

Specifically, it was shown that in d = 2 the estimate D es t(At) of the tracer diffusion constants 
depend on the time interval At of estimation as 

D est (At) = A + B\n(At)/At + C / At , (11) 

where A, B, C are phenomenological constants. Therefore we have analyzed D es t{At) as a function 
of At in the present case (Fig. right part). We found rather generally that there is a significant 
dependence of D es t{At) on At for At < 2.10 3 , while for At > 5.10 3 the dependence on At can 
safely be neglected. A remarkable feature of the results also is that the faster B particles exhibit 
(in the example shown in Fig. [SJ) a diffusion constant that is only about a factor of three larger than 
the slower A particles, while the jump rate is a factor of 100 larger. This fact already indicates 
that there is no straightforward relation between the tracer diffusion constants and the jump rates. 

In the description of collective diffusion, the Onsager coefficients Aaa, Aab, A^ and Abb play 
a central role, since they appear as coefficients in the linear relations between particle currents 
3A>3b and the corresponding driving forces, the gradients of the potential differences between A 
(or B) particles and vacancies V, respectively [171 ]: 

j A = -(A A A/k B T)V(fi A - Mv) " (AAB/k B T)V(vB - Hv), (12) 



j B = -{A B A/k B T)V(nA - Vv) ~ (A B B/k B T)V(nB - Hv) 



(13) 



Note that due to the symmetry relation 



A BA = A A 



B 



(14) 



only three of these four Onsager coefficients are thought to be independent. There is no simple 
relation between the two jump rates Ta,^b (and temperature T and the concentrations ca, cb) and 
these three Onsager coefficients Aaa, A-ab, A-BB, of course. Hence it is a task of the simulation to 



estimate these Onsager coefficients, and it is well known 



17 



that this can be done by applying 
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a force to the particles, which acts in the same way as a chemical potential gradient. Due to the 
periodic boundary conditions, particles that leave the box at one side will reenter at the opposite 
one, and hence a chemical potential gradient causes a steady state flux of particles through the 
simulation box in the direction of this driving force. Care is needed in two respects: 

• One must average long enough to make sure that slow transients after the imposition of the 
force have died out and steady-state conditions are actually reached. 

• One must make sure that the applied force is small enough so one works in the region where 
the response of the system to this force is strictly linear, as written in Eqs. (|12|). (I13|) . and 
nonlinear corrections can be completely neglected. 



This method of estimating Onsager coefficients was pioneered by Murch and Thorn |23] for 



one-component lattice gases and extended to random alloy models in 



17 



3 



. We refer the reader 



Id tlim« papers lor a more del. ailed justitication and discussion of this method, hollowing Q we 
implement this force on species of particles 7 (7 = A or B) by taking the jump rates in the in- 
direction as 

rw = &r 7 , t { 2 ] x = b' 1 r 7 , b > 1 (15) 

while the jump rate in the ±y— directions remains T 7 . If we would have a single particle (s.p.) 
only, the mean velocity in the +x— direction would be v x ' p ' = r 7 (6 — which should correspond 
to v x p ' = (Try/kBT)F x , F x being the force in x— direction, in the regime of linear response. Hence 
one concludes that from the velocity of species a one can deduce the Onsager coefficient A a7 if a 
force F x is exerted on species 7 via 

4 a) K p - = (T a c a )~ l A a ^. (16) 

The application of this method is illustrated in Fig.El There the mean displacement (x) of A— and 
B— particles is followed over 2.5 x 10 4 Monte Carlo steps (MCS) per site, and a very good linearity 
of (x) vs. t is observed (left part). In order to check for nonlinear effects, the bias parameter b is 
varied in the range 1.05 < b < 1.5, and the results are extrapolated to b — ► 1. (right part of Fig. EJ). 
Consistent with previous work on the random ABV model nonlinear effects are rather weak, 
and in this way we are able to estimate Onsager coefficients with a relative error of a few percent. 

Still a different approach was followed to estimate the inter diffusion constant Di n t- We prepare 
a system in thermal equilibrium in the presence of a wavevector-dependent chemical potential 
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difference 6(jl(x) defined as 

2n 

5/j,(x) = ha(%) - Vb(x) = <5cos(— x), (17) 

S being an amplitude that needs to be chosen such that the resulting concentration variation is 
still in the regime where linear response holds, and A is the wavelength of the modulation (which 
is chosen such that the linear dimension L of the lattice is an integer multiple of A). Note that in 
the Ising spin representation 5^x{x) simply translates in a wavelength-dependent magnetic field, of 
course. The system then is equilibrated in the presence of this perturbation for a large number of 
Monte Carlo steps (of the order of 10 6 MCS). This causes a corresponding periodic concentration 
variation, see Fig. [7[ left part. The sinusoidal shape of this initial concentration variation provides 
a confirmation that the linear response description is applicable otherwise the presence of higher 
harmonics in the concentration variation would indicate the presence of nonlinear effects. Then 
a "clock" is set to time t = and the perturbation 5fi(x) is put to zero for times t > 0. As a 
consequence, the concentration variation decays to zero as the time t — > oo. It turns out that 
this decay with time can be described by a superposition of two simple exponential decays, one 
governing the decay of the concentration difference 5c(x) = ca{x) — cb{x) of the particles, the other 
corresponding to the decay of the total density. As discussed in detail for the random ABV model 
[~Lt1 |. the concentration variation can be described therefore as (k = 2ir/\, and D + > D_ are two 
diffusion constants) 

Sc A (t) = c\ exp(-D + k 2 t) + c A exp(-Z)_fc 2 f) , (18) 



5c B (t) = c\ exp(-D + k 2 t) + c~ B exp(-D_k 2 t) , (19) 

where c\, c7, cj, are amplitude prefactors, which one can estimate from the treatment that will 
be outlined in the following section. Here we only mention that c\+c^ = 5ca(0), c^+c^ = 5cb(0), 
and in the limit cy — > we have c\, oc cy — ► 0, while c7, stay finite (of the order of 8). In this 
limit the two diffusion constants D + , D_ are of very different order of magnitude, since -D_ oc cy, 
while D + stays of order unity ^t|. Thus density variations have a very small amplitude (of order 
cy) and decay fast, while concentration variations decay much slower. This consideration leads us 
to identify D_ as the interdiffusion constant Di n t in this limit. For finite nonzero cy, however, 
in principle both density and concentration variations are coupled, and both diffusion constants 
D + , D- contribute to the interdiffusion of A and B particles [13] . 
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The right part of Fig. [7| illustrates that even for cy as large as cy = 0.04 there is already 
a reasonable separation between density and concentration fluctuations: both 5c Ait) and (5ce(i) 
reach their asymptotic decay (where only the same factor exp(— D_k 2 t) matters, as is evident from 
the fact that there are two parallel straight lines on the semilog plot) already at a time t ~ 2000, 
long before the concentration variations have decayed to zero. 



IV. THEORETICAL PREDICTIONS 



A basic ingredient of all analytical theories are the conservation laws for the numbers of A and 
B particles, which lead to continuity equations for the local concentrations CA[r,t),CB{r,t) 

^M + V-L(r,t) = 0, (20) 



+ V.J,(f,,) = 0, (2!) 

Note that these equations hold rigorously, if a local concentration field c a (r, t) {a = A, B} can be 
defined, unlike the so-called constitutive relations, Eqs. p2j l .(|lii| ) . which are only approximately 
true: these equations only are supposed to hold in the case that the gradients VifiA — fJ>v), V(/ie — 
fj,y) are sufficiently small, otherwise the relation between currents and gradients is non linear. In 
addition, a second requirement is that statistical fluctuations can be neglected; otherwise a random 
force term needs to be added on the right hand side of Eqs. (fT2|) . (fT3 |) j^ . We also note that in 
our model (unlike real alloys, where vacancies can be created by hopping of atoms from lattice 
sites to interstitial sites, and where vacancies can be destroyed by hopping of interstitial atoms to 
a neighboring vacant site of the lattice □ □ □ ) also vacancies are conserved, and hence 

^M + V.iv-(?, ( )=0. m 

However, as discussed in jl^l there is no need to include cy(r, t) and jy(r, t) as additional dynamical 
variables in the problem: the condition that every lattice site is either occupied by an A-atom, 
B-atom or vacant (V) translates into the constraint c^(r, t) + ce(r, t) + cy(r, t) = 1. Similarly, one 
finds that jy = ~ij A + 3b) Q- 

In order to be able to relate the chemical potentials in Eqs. ()12j) . ()lH|) to the concentration 
variables, we use the thermodynamic relation 

Ma = (-Of) , (23) 



dN a/ T ,N^ a) 
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N a being the number of particles of species a, and F being the total free energy of the system. We 
decompose F into the internal energy U and the entropic contribution —TS, with S being simply 
the entropy of mixing 

S = -k B [N A In N A + N B In N B + N v In Ny - N In N] , (24) 

where N = Na + N B + Ny then is the total number of sites on the lattice, and c a = N a /N then 
is the concentration of species a. While Eq. 1)24(1 is exact in the non-interacting ABV model, it 
still holds in the disordered phase of the interacting model in the framework of the Bragg- Williams 
mean field approximation. In the disordered phase, no sublattices need to be introduced, and then 
the concentration variables on average are the same for all lattice sites. Then U can be written as 

U = -Nz{eAAc\ + 2e AB c A c B + e BB c 2 B ), (25) 

where z is the coordination number of the lattice, and consistent with the simulated model (Sec. II) 
a nearest neighbor interaction is assumed. Note that the basic approximation of Eq. (|25[) is the 
neglect of any correlation in the occupancy of neighboring lattice sites. 



With some algebra 
equations 



13] one can reduce Eqs. ([12)1 . (fl3|) . ([20 )1 - (|25 |) to a set of two coupled diffusion 



^ = J> Q/3 V%, (26) 



where the elements D a p of the diffusion matrix are given by 



Daa = Aaa(- + - + ^) + Aab(- + ^), (27) 

CA Cy k B T Cy k B T 



Dab = Aaa(- + ) + A AB (- + - + ^) , (28) 

Cy k B l C B Cy k B l 



Dba = Aab(~ + - + f^) + Abb(- + ff ) , (29) 

CA Cy k B l Cy k B l 



D BB = A AB (±- + ) + A BB (- + - + ) ■ (30) 
cy k B T c B cy k B T 

Note that Da B ^ D B a- Introducing Fourier transforms and diagonalizing the diffusion matrix the 
solution indeed can be cast into the form of Eqs. (|18|) . (|19|) . As has already been mentioned in this 
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context, for cy — > the two eigenvalues D+,D- of the diffusion matrix adopt very different orders 
of magnitude jl 71 ] : 

D+ « (A AA + 2A J 4£ + A BB )/c V , (31) 



g _ ^ A AA A BB - 1 | 1 2ze 

Aaa + 2Aab + A B b ca cb k B T 

Since in this limit the A Q( g oc cy, the coefficient D + reaches a finite limit for cy — > 0, while 

-D- oc cy- We also recognize that D_ can be decomposed into a product of two factors: a "kinetic 

factor" Ai n t, composed by a combination of Onsager coefficients, and a "thermodynamic factor", 

which is nothing but an effective inverse "susceptibility" x~ l describing concentration fluctuations, 

normalized per lattice site, 



X 



1 = c A x + c B l - 2ze/k B T = [cx(l - caT 1 - 2ze/k B T (33) 



In the last step, we used the fact that eg = 1 — ca for cy — > 0. We call % a "susceptibility" because 
in the translation to the Ising spin representation \ simply becomes proportional to the derivative 
of the "magnetization" with respect to the field. Note that for e > (i.e., a mixture with unmixing 
tendency) Eq. (|33|) exhibits a vanishing of x an d hence of the interdiffusion constant D_ at the 
mean field spinodal curve, defined by 

k B T s (c A )/e = 2cx(l - c A )z . (34) 

The mean field spinodal touches the coexistence curve of such a phase-separating mixture at its 
maximum in the critical temperature, i.e. ksT^ F /e = z/2 = 2, for a square lattice. Actually, the 
symmetry of the Ising Hamiltonian in zero field implies that the maximum critical temperature of 
the Ising antiferromagnet, which occurs at zero field as well, then is also given by 

k B T c % x /\e\ = z/2 = 2, e<0 (35) 

Comparing this estimate to the exact result, Eq.(|2J), we notice that the mean field approximation 
actually overestimates the maximum critical temperature of the ordering alloy by almost a factor 
of two, as is well known. Note that this error increases for ca ^ 1/2 |37|. So Eq. (|32[) cannot be 
assumed to be quantitatively reliable. Note that for ordering alloys (where e < 0) the interdiffusion 
constants gets enhanced (rather than reduced, as it happens for alloys with unmixing tendency) as 
an effect of the interactions. Beside that, Eq. (|32|) does not predict any singularity of D_ as one 
approaches the order-disorder phase boundary T c (ca) from the disordered side. 
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Discussing now the kinetic factor Aj n t, we recall the popular approximation to neglect the 
off-diagonal Onsager coefficient in comparison to the diagonal ones. This leads to 

A int «A^ + A^ (36) 

With this approximation, Eq. (|32j) reduces to the well-known "slow mode theory " of interdiffusion, 
which has been much debated in the case of fluid polymer mixtures j^-j^. A mean field type 



approximation for self-diffusion 



17, 



42 



43, 



44 



45l | then relates the Onsager coefficients A^, Abb 



and tracer diffusion coefficients Df,Df, 

A AA = Dfc A , Abb = D?c B (37) 

and thus the "slow mode" theory predicts the following relation between tracer diffusion coefficients 
and the interdiffusion constant (remember eg = 1 — ca for cy —* 0) 

D s in T = {(DfcA)- 1 + [Df{\ - c A )]- l }{[c A (l ~ ca)]- 1 - 2ze/k B T}. (38) 



A rather different result, the so-called "fast mode" theory 



46, 



471 ] . can be obtained by several 



distinct arguments. We mention only one of these arguments here, which starts from the assumption 
that everywhere the vacancy concentration cy(f,t) is in thermal equilibrium, i.e. 

V/x„ = 0. (39) 

Of course, in our model Eq. (|39|) cannot be justified, in view of the constraints cy(f, t) = 1 — 
ca(t, t) — cs(r, t), jy = —{ja + 3b) and Eqs. (f2^]) - (|2^|) there is no freedom to make additional 
assumptions on [iy at all, V /J,y(r,t) already is determined from these other equations. However, 
the motivation for Eq. (|39|) is that for real systems there is no strict conservation for the number of 
vacancies: in real (three-dimensional) alloys, vacancies can be created and destroyed by formation 
or annihilation of interstitial atoms, or by interaction with other lattice imperfections such as 
dislocations, grain boundaries, etc. For two-dimensional surface alloys |2]|, vacancies can be created 
and destroyed if an atom from the considered surface monolayer becomes an adatom on top of 
this monolayer, or an adatom executing surface diffusion [18L 1.3 1| becomes incorporated into the 
monolayer via a jump to a vacant site inside the monolayer. In view of these physical mechanisms 
which are forbidden in our model, Eq. ()39|) may represent a physically interesting limiting case. A 
priori, it is not clear for a particular system, whether for the time scales of interest it is closer to 
a situation where vacancies are in equilibirum (Eq. (|39|0 or conserved (Eq. (|22j0 . Our numerical 
studies are concerned with the latter case exclusively. Nevertheless it is of interest to mention that 
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Eq. (|3~9*)) yields also a structure D_ = h-intX 1 but with Ai nt = csDf + CADf and hence one finds 
instead of Eq. (3$ [12] 

S£T = K 1 - c A)Df + c A D?){[c A {l - ca)]- 1 - 2ze/k B T}. (40) 

While for Df » L> t A (a case expected if V A C T^, as used in our simulation) one expects from 
Eq. (|4()|) that the faster diffusing B species dominates interdiffusion, the opposite is true according 
to Eq. IJ38JI : therefore the names "fast mode" and "slow mode" theory have been chosen. In 
both equations (and in Eq. (|32|) . where the off-diagonal Onsager coefficient is not neglected, unlike 
in both these theories) the thermodynamic factor is treated by a simple Bragg- Williams mean 
field approximation, however, which is no problem for the random alloy ABV problem treated in 
Ref. jl7|, but clearly will introduce additional shortcomings in the present case. 

As a final disclaimer of this section we emphasize that Eqs. (|20j) - (|4U[) were meant to provide a 
brief review of "chemical diffusion" (or "collective diffusion" ) in the context of the present lattice 
gas model only, and hence many interesting and important facets of this topic have not been 
mentioned at all and we direct the interested reader to the rich literature on this subject ODD 

□□□□□□□□ 



V. SIMULATION RESULTS 
A. Tracer diffusion 



We start with a discussion of the tracer diffusion coefficients (Figs. |H1 The simplest case 
refers to equal jump rates E4 = Tb = 1 of both types of particles A and B (Fig. El). In the infinite 
temperature limit then there is no longer any physical difference between A and B particles, they 
simply differ only by their labels: then Df = Df , and become independent of concentration ca 
(thick horizontal straight line in Fig. |HJ). Note that for ca = 0.96 there are no B particles since 
cy = 0.04 and then becomes independent of temperature, similarly as Df becomes independent 
of temperature for ca = 0. Of course, the curves for Df are simply the mirror images of those for 
Df around the symmetry line <?Amaix = 0- ~ c v)/2 = 0.48 of the static phase diagram, Fig. |2h>, 
since an interchange of A and B means that ca gets replaced by 1 — cy — cb- 

It is seen that the onset of ordering depresses self-diffusion very strongly, while short range 
order (as it occurs for T = 1.2) has a minor effect only. For T = 0.6, however, the ordering near 
ca = 0.48 is rather perfect and there deep minima of Df, Df occur, the tracer diffusion coefficients 
decrease by about two orders of magnitude. Of course, since Df, Df are not symmetric around 
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c Amax = c B l max = 0- ~ c v)/2 = 0.48, due to the choice of a kinetic Monte Carlo algorithm 
which lacks the symmetry between the motion of an A particle, mediated by a vacancy, in a B 
environment and in an A environment at finite temperatures, the minimum of Df does not occur 
precisely at cj^^,, as is seen from Fig. |H] (left part). In our algorithm, an A particle jumps to a 
vacant site with a jump rate exp(— An\eAB\/kBT) when the difference between the number of 
AB bonds involving an energy €ab each between the initial and final state is An > 0, and with 
a jump rate Fa else. It is easy to be convinced that this algorithm satisfies detailed balance with 
the canonic equilibrium distribution, as it should be ^^]. In the limit ca — > 1 — ey, we always 
have An = 0, so there is no temperature dependence. In the limit ca — ► 0, however, every A 
atom not having a vacancy as nearest neighbors will have four B neighbors on the square lattice, 
while an A atom with a vacancy neighbor has only three B neighbors. As a result, the jump of 
an A atom that has a B neighbor, to a vacant site involves "breaking" an AB bond, and hence 
this rate is suppressed by a factor exp(— {eAsl/hsT). This effect is responsible for the temperature 
dependence of Djr for ca — * 0. 

Kehr et al. 17] presented arguments to relate the tracer diffusion coefficients to Onsager 
coefficients which take a simple form in the case of identical jump rates r^,r^, namely 

= Aaa/ca- ^ba/cb, D? = A bb /cb - Aab/ca (41) 

Using our estimates for the Onsager coefficients at T = 0.6 (see below) in Eq. (|41j) . one sees that 
the trend of the concentration dependence of is reproduced rather well. However, one should 
note that the derivation of Eq. (|41|) is rigorous only for the special case eAB/ksT = 0, because only 
then the distinction between A and B particles forming the environment of a tagged A particle 
can be neglected. 

When r,4 ^ Tb the self-diffusion coefficients Df and Df lack any symmetric relation of their 



concentration dependence already in the random alloy limit [17J, and for eab/^bT < we are not 
aware of any theoretical treatment to which our simulation results (Fig. E3) could be compared. 
Interestingly, for not too low temperatures (such as T = 0.912, T = 1.2) the concentration de- 
pendence of (the slower diffusing species, since Ta/^b = 0.01 has been chosen in Fig. EI) is 
rather weak throughout, while for Df we have a strong decrease when ca increases up to about 
C A max = 0.48. For ca > c C Amax a g a i n a very weak concentration dependence results. For T = 0.6 
again pronounced minima near c^"!l are seen. Now, for Df we have a strong decrease when ca 
increases up to about c^pj^ aa; , while for ca > c C Arnax a S a i n a vei T weak concentration dependence 
results. 
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Moreover, when for ca = c C Amax ^ ne or der of the AB— checkerboard structure is perfect (apart 
from a four per cent of vacant sites in the system), a jump of an atom to a vacant site occurs with 
rates exp(3eABAB^) or Tb exp(3eAB/^BT), respectively, while the backward jump occurs at 
rates Ta,^b- As a result, a high probability for backward jumps is expected, and this is borne out 
by a study of the correlation factor / for self-diffusion (Fig. 1101 right part). Following standard 



treatments 



we decompose tracer diffusion coefficients Dt as 

A = VWf , (42) 



where V is the vacancy availability factor already defined in Eq. Q, and W is the average jump 
rate for the considered particle species. W is easily estimated in the simulation from the ratio 
of the number of performed jumps to the number of all attempted jumps. The product VaWa 
is plotted in Fig. ^] (left part) versus ca at various temperatures. For ^abI^bT — > we simply 
expect a horizontal straight line, VaWa = 0.04, since then Wa = 1, Va = cy (a\ = 0). There 
is no independent way to determine /, however. Therefore Eq. (|42jl is taken as a definition of /, 
to be derived from Dt, while the tracer diffusion constants are estimated from the mean square 
displacements of the tagged particles, as explained in Sec. IIIII of this paper. For ca — > 0.96, 
when no B particles are present, the temperature dependence drops out and / reduces to the 
value / = 0.487 known from studies of a one-component non-interacting lattice gas on a square 
lattice with concentration ca = 0.96 29]. Note that our data for Dt, V, W and / at the hi gher 
temperatures (where no order-disorder transition occurs) resemble analogous results of Murch [26j 
for a simple cubic alloy. 

As a final comment about self-diffusion, we consider the temperature dependence of Df and 
Df for the critical concentration c^ 4 * = c c ^ lt = 0.48 (Fig. 111(1 . One sees that at high temperatures 
(T > 2T C ) the temperature dependence is very weak, and the tracer diffusion coefficients settle 
down at their infinite temperature asymptotes. Approaching the critical point one sees a more 
rapid decrease of both Df and Df , with a maximum slope presumably right at T c , while for T 
below T c a crossover to the expected thermally activated behavior at low temperatures occurs. In 
fact, one expects that Dt — D* oc (T — T C ) 1_Q [2^, where D* is the value of the tracer diffusion 
coefficient at the critical point, and a is the specific heat exponent of the model. However, for the 
two-dimensional Ising model a = i.e. the specific heat has a logarithmic singularity 

only. The insert of Fig. 1111 shows a log-log plot of Dt — D* versus (T — T c )/T c , and one sees that 
the data are compatible with a power law with slope of unity; presumably the accuracy of our 
simulations does not suffice to identify the presence of a logarithmic singularity in our data. 
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B. Onsager coefficients 

As a first issue of this subsection, we turn to the concentration dependence of the Onsager 
coefficients (Figs. El ECU)- For Ta = ^b ah Onsager coefficients are symmetric around ca = cb = 
(1 — cy)/2, as it must be, while for Ta they are not. We have also included an approximate 

relation suggested by Kehr et al. [lj| between Onsager coefficients and tracer diffusion coefficients, 
namely 

l-/(c) cpDl 



A a /3 = c a Dt 



$af3 + 



/(C) ECyA 1 



(43) 



where c = ca + cb, f(c) being the correlation factor for tagged-particle diffusion in a lattice gas 
with summary concentration c. It is seen that this relation accounts for the general trend of the 
diagonal Onsager coefficients rather well, although for the off-diagonal Onsager coefficient it seems 
to work only qualitatively (Fig. I13JI . In the regime of the ordered phase the diagonal Onsager 
coefficients (note the logarithmic ordinate scale) are distinctly smaller than for ca — ► or cb — > 0, 
respectively, when Fa = ^b- 

An interesting aspect of the off-diagonal Onsager coefficient Aab = Aba (Fig. Il3|) is that it is 
essentially zero for ca — ► if Ta = ^b while for Ta/^b = 0.01 it is essentially negative in this 
limit. A negative Onsager coefficient means that the currents of A and B particles are oriented 
in the opposite direction. A further change of sign of this off-diagonal coefficient is found near 
the phase boundary of the order-disorder transition; but near ca = cb = (1 — cy)/2 the Onsager 
coefficient seems to be positive again, although its absolute value seems to be very small. We do 
not have any clear physical interpretation for this surprising behavior. Note also, that Eq. I|43|) can 
never yield a negative Onsager coefficient, since < /(c) < 1 by definition, and hence all terms in 
Eq. 03] are non- negative. 

Finally Fig. El shows the temperature dependence of the Onsager coefficients for the concen- 
tration ca = cb = (1 — cy)/2 where the critical temperature T c of the order-disorder transition is 
maximal. Note that for Ta/^b = 0.01 the magnitude of the off-diagonal Onsager coefficient Aab 
is comparable to the smaller (Aaa) of the diagonal ones, both at very high and at very low tem- 
peratures. This finding confirms the conclusion of Kehr et al. [121], that in general the off-diagonal 
Onsager coefficient must not be neglected. We also note that the general trend of the temperature 
dependence of the Onsager coefficients is very similar to the behavior of the self-diffusion coeffi- 
cient, see Fig. Both quantities reflect the strong decrease of mobility of the particles at low 
temperatures. 



19 



C. 



Interdiffusion 



Fig ED presents a plot of the interdiffusion constant Dj„[ vs. concentration for the case of 
equal jump rates (r^ = r# = r = l)atT = 0.6 and compares the results to various analytical 
approximations: D_ (Eq. (J22J)), the "slow mode" expression Df n ™- (Eq. (JHSJ)), the "fast mode" 
expression D^' (Eq. (|4Uj)). and a very simple result justified by Kehr et al. jl^] for the non- 
interacting random alloy model, 



While this last expression overestimates the numerical results, all other expressions underestimate 
them significantly. It is seen that in this case there is not much difference between the slow mode 
and fast mode theory, but both are off from the data. In this case using the full expression (Eq. (|32[l) 
presents no improvement, unlike the non-interacting case. Of course, at finite temperature in d = 2 
the mean field theory implicit in Eq. (|32|) is not expected to be accurate at all. 

It now is no surprise any longer that in the asymmetric case Ta/Tb = 0.01 the various ap- 
proximate expressions are not reliable either (Fig. Ilfi|h In particular, for concentrations near 
ca = cb = (1 — cy)/2 = 0.48 a pronounced minimum is predicted, while the actual simulation re- 
sults reveal a rather shallow minimum only. Again the conclusion is that there is no reliable simple 
relation between self-diffusion and interdiffusion coefficients, and the temperature dependence of 
Di n t at ca = cb = 0.48 at higher temperatures (Figs. 117118)) confirms this conclusion. Again, for 
Fa = r# = T = 1 Eq. (|H)) is closest to the data, while Eq. (|52~)) is worst. For T — > oo, however, 
in this limit for ca = cb = (1 — cy)/2 and Ta = r# = 1 all expressions coincide (at a value 
highlighted by an arrow in Fig. 117)1 . and the numerical data have been found in good agreement 
with this prediction 17]. Thus it is clear that including interactions among the particles destroys 
the applicability of the simple theories. 



In this paper, the study of mobility of particles, interdiffusion and tracer diffusion coefficients 
of a lattice model for a binary alloy, that was presented in Ref. for the simple non-interacting 
limit only, has been extended to the case where an attractive nearest-neighbor interaction between 
unlike particles leads to an order-disorder transition on the considered square lattice. While most 
theoretical considerations of the previous work can be simply extended to the present case, the 
mean-field character of the approximations that are involved clearly emerges as a severe limitation 




(44) 



VI. CONCLUSIONS 
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of the usefulness of all these approaches. On the other hand the Monte Carlo techniques described 
in Ref. suitable for the direct estimation of all Onsager coefficients and the interdiffusion 
constant Di n t as function of the ratio of jump rates Ta/Tb, temperature T and concentration ca, 
are rather straightforward to apply. Exploring this rather large parameter space numerically is 
however somewhat tedious, and an understanding of diffusion phenomena within the framework 
of lattice models for interacting particles by simple analytical expressions clearly would be desir- 
able. However, the approximate expressions discussed in the present paper clearly do not give 
qualitatively accurate results. 

Of course, the present study is a first step only: in order to make closer contact with possible 
experiments in surface layers of metallic alloys, it would be interesting to consider other lattice sym- 
metries (triangular and centered rectangular lattice rather than square lattices), further neighbors 
interactions, etc.. 

A very important extension would also be the inclusion of asymmetric effects (caa € bb) and 
nonzero energy parameters involving vacancies (€av,€ rv, evv)- Thus effects could be described 
that vacancies occupy preferentially sites at interfaces \l\ or in one of the sublattices (2^ . Such 
effects are expected to modify the diffusion behavior significantly. 

We thus hope the present study will stimulate the development of more accurate theoretical 
descriptions of diffusion phenomena in alloys that undergo order-disorder transitions. Also corre- 
sponding experiments studying a wide range of temperature and composition, would be desirable. 
Then it might be worthwhile to combine the present kinetic Monte Carlo methodology with "ab 
initio" calculation of jump rates, ordering energies e a p, etc. 
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FIG. 1: Schematic view of the (100) surface of a substrate (shown as large open circles), whose periodic 
potential provides a square lattice of prefered adsorption sites (which here are assumed in the center of the 
square formed by the substrate atoms), yl-atoms are shown as black circles, -B-atoms are shown as grey 
circles, and vacancies (V) are shown as empty circles. The energies of the nearest-neighbor interactions 
between different kind of atoms (indicated by thick lines) are labeled by caAi £bb and cab respectively. 
The simple choice 6ab = e i £ aa = e BB = is taken througout. This means that A— atoms prefer B— atoms 
as nearest neighbors, but it does not matter whether its nearest neighbors are also A— atoms or vacancies, 
respectively. The jump rates for A — V and B — V exchanges are labeled by Ta and Tb, respectively. For 
simplicity, the B— atoms are considered as the faster particles (r^ = 1), and Ta < Tb- 
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FIG. 2: (a) Phase boundary for the order-disorder transition of the ABV model with cv = 0.04. The 
phase boundary of the pure Ising antiferromagnet [3^ (equivalent to the case cv = 0) is also included for 
comparison, as a dashed line, (b) Critical curve T vs. ca, where ca = 1 — cv — cb- The critical curve 
of the pure Ising antiferromagnet is also included for comparison, as a dashed line. Temperature T is 
always measured in units of the maximal critical temperature T c of the pure model (no vacancies, cy = 
and ca = cb = 0.5), cf. text. 
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FIG. 3: Dependence of the staggered susceptibility \ (a) and the fourth order cumulant Ul (b) on the 
temperature, along the critical line H = corresponding to the critical concentration ca = 0.48. Several 
system sizes are considered, as indicated. Here T c denotes the maximal critical temperature of the model 
without vacancies (H = then corresponds to ca = 0.5). 
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FIG. 4: (a) Plot of the size-dependent critical temperature T C (L) (denned as the maximum of x(T, L)), in 
terms of the scaled variable l}l v . The critical Ising exponent v = 1 is employed. The linear extrapolation 
to the thermodynamic limit, shown as a dashed line, provides an estimation of T c (c,y = 0.04) = 0.905(5) 
for the ABV model with cy = 0.04. (b) Scaling plot of the order parameter ip. The estimated critical 
temperature and the Ising critical exponents v = 1 and (3=1/8 are employed. 
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FIG. 5: Determination of the tracer diffusion coefficients, (a) Mean square displacements of tagged A 
particles (open dots) and B particles (full dots) as a function of Monte Carlo steps per particle. The 
temperature is T — 1.2 (in units of the Ising critical temperature), and the concentrations are ca = 0.4, 
cb = 0.56. The ratio of jump rates is Ta/^b = 0.01. (b) Estimates of the tracer diffusion coefficients as a 
function of the time interval used. The lines represent the fits of the data after using Eq. Qlip. 
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FIG. 6: Determination of the Onsager coefficients, (a) Mean displacements along the x— direction of A 
and B particles as a function of Monte Carlo steps per particle. The temperature is T = 0.6, and the 
concentrations are ca = 0.71, cb = 0.25. The ratio of jump rates is Ta/Tb = 0.01, and the bias parameter 
is b = 1.1. (6) Estimates of the Onsager coefficients Ay/cj by extrapolation to bias parameter 6 = 1. 
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FIG. 7: Determination of the interdiffusion coefficients. For t < we impose a cosine-like varying bulk 
field H(x) which introduces a modulation in the concentration of A and B particles. The characteristic 
length of this perturbation is A = 32 lattice spacings. (a) Temporal evolution of the concentration of A 
particles along the x— direction in the lattice. Times correspond to t = (circles), t — 10000 (squares) 
and t — 20000 (diamonds), respectively. The thick line marks the wavelenght A = 32 of the applied bulk 
field. The temperature is T = 1.5, and the concentrations are ca = Cb = 0.48. The ratio of jump rates 
is Ta/^b = 0.01. (b) Amplitude of concentration profiles as a function of time, for A particles (circles) 
and B particles (squares) . The dashed lines correspond to fits of the data to single exponential functions, 
characterized by a decay constant D int . 
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FIG. 8: Tracer diffusion coefficients for A particles (left) and B particles (right), as a function of the 
concentration ca- The jump rates are Ta = = 1 and several temperatures are considered: T = 0.6 
(dotted line), T = 0.912 (dashed line) and T = 1.2 (dot-dashed line). The thick line indicates in both cases 
the noninteracting, infinite temperature limit (random alloy model). Dots represent results obtained from 
Eq. gTJ. 
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FIG. 9: Tracer diffusion coefficients for A particles (left) and B particles (right), as a function of the 
concentration ca- The jump rates are Ta/Tb = 0.01 and several temperatures are considered: T = 0.6 
(dotted line), T = 0.912 (dashed line) and T = 1.2 (dot-dashed line). The thick line indicates in both cases 
the noninteracting, infinite temperature limit (random alloy model). 
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FIG. 10: Effective jump rate (left) and correlation factor (right) for A particles, as a function of the 
concentration ca and for different temperatures: T — 0.6 (dot line), T — 0.91 (dashed line) and T = 1.2 
(dot-dashed line). The jump rates are Ta = = 1- 

The first quantity provides an idea of the rate at which jumps actually occur at a certain temperature and 
composition. It is defined as the product of the vacancy availability factor Va (related to the short-range 
order parameter on) and the average jump rate W a (defined as the quotient of the number of performed 
jumps to the number of all attempted jumps). The T — > oo limit is given by V = c v = 0.04, because in this 
case oti = 0. 

Once we obtain VW we can estimate the correlation factor / applying the definition D t = VWf and using 
D t from FiglHl See Refs. [3 13 13 13 for details on the effect of correlations on tracer diffusion in lattice 
gas models. 

The limit value / = 0.487 for ca — > 0.96 is known from Ref. 29j. This corresponds to a noninteracting, 
one-component lattice gas in a square lattice with concentration c = 0.96. 
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FIG. 11: Dependence of the tracer diffusion coefficients on the temperature, for a stoichiometric composition 
ca = cb = 0.48. The ratio of jump rates is Ta/^b = 0.01. Circles are A particles and squares are B particles. 
The dashed arrow marks the critical temperature T c = 0.905 (in units of the Ising critical temperature), 
while the thick arrows indicate the asymptotic, infinite temperature values for both coefficients. 
Inset (up): Arrhcnius plot of D t for T < T c . Inset (bottom): scaling plot of \D t - D*\ ~ |T - T^ 01 with 
a = (specific heat exponent of the Ising model) . The dashed line has a slope of unity. 




FIG. 12: Plot of the Onsager coefficients A« as a function of the concentration ca, for a fixed temperature 
T = 0.6. The jump rates are Ta = = 1 (left) and Ta/^b = 0.01 (right). The lines correspond to data 
obtained directly from the simulations, while the points correspond to the estimates obtained after using 
the aproximation of Eq. (|43[) for the random alloy model. 
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FIG. 13: The same as figureElbut here the crossed coefficient Ay is shown. 
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FIG. 14: Plot of the Onsagcr coefficients Ajj as a function of the temperature, for a stoichiometric concen- 
tration ca = cb = 0.48 and jump rates Ta/^b = 0.01. The lines are drawn to guide the eyes. 
Inset: plot of the same coefficients in the vicinity of the critical temperature, showing the linear approach 
to finite values right at T c . 
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FIG. 15: Plots of the intcrdiffusion coefficient as a function of the concentration ca, for a temperature 
T = 0.6. The jump rates are Ta = = 1- The arrows mark the corresponding critical values of ca for 
this temperature. Simulation results are compared to different theoretical approaches, for a discussion see 
the text. 




FIG. 16: Plots of the intcrdiffusion coefficient as a function of the concentration ca, for a temperature 
T = 0.6. The jump rates are T A /T B = 0.01. 




FIG. 17: Logarithmic plot of the interdiffusion coefficient as a function of the temperature, for a concentra- 
tion ca = Cb = 0.48. The jump rates arc Ta = ^b = 1- The arrow marks the infinite temperature result, 
where all the quantities showed in the plot coincide. 




FIG. 18: Logarithmic plot of the interdiffusion coefficient as a function of the temperature, for a concen- 
tration ca = cb = 0.48. The jump rates are Ta/^b = 0.01. 



